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Abstract 

Consider a model of particles (nucleons) which has a two-body interaction 
which leads to bound composites with saturation properties. These properties 
are : all composites have the same density and the ground state energies of 
composites with k nucleons are given by —kW + ak"^/^ where W and a are 
positive constants. W represents a volume term and a a surface tension term. 
These values are taken from nuclear physics. We show that in the large 
limit where is the number of particles such an assembly in a large enclosure 
at finite temperature shows properties of liquid-gas phase transition. We do 
not use the two-body interaction but the gross properties of the composites 
only. We show that (a) the p — p isotherms show a region where pressure 
does not change as p changes just as in Maxwell construction of a Van der 
Waals gas, (b) in this region the chemical potential does not change and (c) 
the model obeys the celebrated Clausius-Clapeyron relations. A scaling law 
for the yields of composites emerges. 

For a finite number of particles A^ (upto some thousands) the problem can 
be easily solved on a computer. This allows us to study finite particle number 
effects which modify phase transition effects. 

The model is calculationally simple. Monte-Carlo simulations are not 
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needed. 
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I. INTRODUCTION 



A very popular and highly successful model for collisions of two nuclei at intermediate 
energies (50 MeV to 100 MeV per nucleon) is the following. Because of many collisions 
between nucleons, a statistical equilibrium is reached. The temperature rises. The system 
expands from normal density and composites are formed on the way to disassembly. As 
the system reaches between 3 to 6 times the normal volume, the interactions between the 
composites become unimportant (except for the long range Coulomb interaction) and one 
can do a statistical equilibrium calculation to obtain the yields of the composites at a 
volume which is called the freeze-out volume. Although the model is simple, actual reahstic 
calculations based on the model are much harder. The nucleus is a finite system. It has 
two kinds of particles, neutrons and protons (generically termed nucleons). Protons carry 
charges and prevent large nuclei from being formed. For realistic treatment, the idea of a 
strict freeze-out volume has to be modified. 

Here we consider the same physics but with the following simplifications: only one kind 
of particle is considered and the Coulomb interaction is neglected meaning arbitrarily large 
"nuclei" can be formed. The energy scale is MeV (million electron volt) and the length scale 
is fm (lO^^^cm) so the sahent features of nuclear physics are retained. The binding energy 
and the volume of a composite is proportional to the number of particles (nucleons) in the 
composite and have a surface tension proportional to the surface area. We show, with rather 
little effort, that the model leads to a first order phase transition as either the density the 
temperature or both are varied. The system has a region of hquid-gas coexistence where, 
as for the Maxwell construction of a Van der Waals gas, pressure remains constant when 
the density increases along the isotherm. In this region the chemical potential remains 
unchanged. As one traverses the path from the liquid phase to the gas phase the Clausius- 
Clapeyron relationship is obeyed. For large systems, a scaling law for composites emerges: 
if we know the yields of composites for one large system, we know these for another large 
system. 
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A more realistic version of this model has been used for Bevalac physics (>250 MeV per 
nucleon beam energy in the lab) by many authors more than twenty- five years ago. It is not 
possible to quote all the references but a review article [1] has a more complete list. The 
possibility of a phase transition was not considered as the collision energies were too high 
for the liquid phase and only very light composites could be produced. 

Phase transitions in heavy ion collisions at intermediate energies became a topic of 
considerable interest starting from the mid-eighties and continues to be a central issue. 
There are many approaches and a large literature too numerous to list. We will refer here 
to only a few which closely follow the underlying physics considered here. The same model 
as used here was adopted in [2] for finite nuclei. By extrapolation it was shown that the 
model leads to a first order phase transition. A brief application of this model is given 
in section VII. A grand canonical model was used in [3] which demonstrated a first order 
phase transition. The approach was quite different from what is used here. We use simpler, 
more traditional and numerical methods. Our results are similar but sufficiently different to 
warrant a full description. A discussion of Clausius-Clapeyron relations and a scaling law 
highlight some interesting physics. 

The celebrated statistical Multifragmentation model (SMM) of Copenhagen [4] , the mi- 
crocanonical models of Gross and Randrup, Koonin [5,6] use the same underlying physics 
as in this work. But the emphasis was on trying to be as close to to the actual nuclear 
situation as one can and thus the phase transition aspects are largely hidden. 



If we have Ua particles of type a, rib particles of type b, ric particles of type c etc. all 
enclosed in a volume V and interactions between particles can be neglected, the grand 
partition function for this case can be written as 



II. THE BASIC FORMULAE 




n 




(1) 



i=a,b.,c... 
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Here the is the chemical potential and Zi the canonical partition function of one particle 
of type i. The average number of particles of type i is given by d{lnZgr)/d{Pni) : 



It is possible that one of the species can be built from two other species. In reverse, a heavier 
species can also break up into two lighter species. If a number of particles of type a can 
combine with P number of particles of type b to produce 7 number of particles of type c, 
then chemical equilibrium implies [7] that the chemical potentials of a, 6 and c are related 

by + (3iJ,b = 'flJ^c- 

In our model we have N nucleons in a volume V but these nucleons can be singles 
or form bound dimers, trimers etc. Chemical equilibrium implies that a composite with 
k bound nucleons has a chemical potential k/i where is the chemical potential of the 
monomer (nucleon). Thus our ensemble has monomers, dimers, trimers etc. upto some 
species with kmax bound nucleons where ideally kmax oo- For practical calculations, we 
use a finite value of k^ax- Most of the results shown here use k^ax = 2000 although we 
have also done calculations with much larger values. Choosing k^ax — 2000 does not mean 
that the total number of nucleons is 2000. The total number of nucleons can be infinite 
but the largest species allowed in the calculation (is somewhat artificially) limited to 2000. 
The total number of nucleons will be denoted by N where N is very large. The quantity 
kmax plays an essential role, setting k^ax too low (for example 200 as shown in section VII) 
makes then liquid-gas transition disappear. An assembly with ~200 particles or less does 
not display the typical behaviour of liquid-gas co-existence. 

We now look into Zi, the partition function of one composite of i nucleons.This factors into 
two parts, a traditional translation energy part and an intrinsic part: Zi = Zi{tran)zi{int) 
where 



n, 




(2) 




(3) 



V_ 
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The intrinsic part Zi{int) of course contains the key to phase transition. If we regard each 
composite to exist only in a ground state with energy ef^, then Ziiint) = exp{—(3ef^). We 
use ef^ = —iW + ai'^^^ where nuclear physics sets 1^=16 MeV and a = 18 Mev. This simple 
model itself will lead to the main results of this paper. Because of the surface term, energy 
per particle drops as i grows. Let us denote by F the free energy of the nucleons where N 
is the total number of nucleons; E be the energy and S, the entropy: F = E — TS. At finite 
temperature F will go to its minimum value. The key issue is how the system of N nucleons 
breaks up into clusters of different sizes as the temperature changes. At low temperature 
E and hence F minimises by forming very large clusters (liquid). But as the temperature 
increases S will increase by forming larger number of clusters thus breaking up the big 
clusters. Gaseous phase will appear. How exactly this will happen requires calculation and 
we find that the system goes through a first order liquid-gas phase transition. 

We used here a shghtly more sophisticated model for Zi{int). This does not make the 
calculation any harder (or alter the qualitative features) but makes it more realistic. We 
make the surface tension temperature dependent in conformity with usual parametrisation 
[%aiT)^ao[{T^-T')/{T^ + T')]'/' 

Here ctq =18 MeV and Tc— IS MeV. At T = Tc surface tension vanishes and we have 
a fluid only. For us this is unimportant as our focus will be the temperature range 3 to 8 
MeV. 

In computing the partition function Zi{int) we include not just the ground state but also 
excited states of the composite in an approximate fashion. We should compute Zi{int) = 
exp(— /Jef^) + JJ29i{^) exp(— /?e). Here e > ef^ and gi{e) is the density of excited states of 
this particular composite. Instead of trying to calculate Zi{int) by performing the sum and 
integral we use a well-known trick. Utilise the relation Zi{int) = exp{—fi{int)/T) where 
fi{int) — ef — Tsi and now use the Fermi- gas formula for the nucleus with i nucleons 
(approximately correct and widely used at intermediate temperature). This gives ef = 
ef^ + iT'^/eo- This is similar to electron gas at finite temperature (excitation energy goes 
like T^) except that in nuclear physics the value of eo is ~16 MeV. The intrinsic entropy of 
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the nucleus at this temperature is 2iT/eQ. The expression for Zi{int) is now complete and 
easily tractable. 

Let us now summarise the relevant equations. For k = 1 (the nucleon which has no 
excited states) 

ni = ^(27rmr)^/2exp(///T) (4) 

and for A; > 1 

Uk = ^{27rmTf/^k^/^ exp[{iik + Wk + kT^eo - a{T)k^/^) /T] (5) 

Here is the average number of composites with k nucleons. In the rest of the paper, for 
brevity, we will omit the qualifier "average" . It is always implied. 
A useful quantity is the multiplicity defined as 

^max 

M=Y.nk (6) 
fe=i 

The number of nucleons bound in a composite with k nucleons is kuk and obviously N = 
Y!k^i kuk- The pressure is given by 

P=Y.^T (7) 

k=l ^ 

Quantities like V, Uk are all extensive variables. These equations can all be cast in terms 
of intensive variables hke N/V — p,nk/N etc so that we can assume both N and V approach 
very large values and fluctuations in the number of particles can be ignored. Thus for a given 
temperature and density we solve for using 

p = (2^^p'^' (exp(/x/r) + 'f]' A;5/2exp[(/xA; + Wk + kT^e^ - a{T)k^/')/T]) (8) 

fc=2 

The sum rule N — Y^lzT ^^fe changes to 1 — J2 kuk/N. 

Prom what we have described so far it would appear that V in eqs.(3), (4) and (6) is 
the freeze-out volume V, the volume to which the system has expanded. Actually if the 
freeze-out volume is V then in these equations we use V which is close to V but less. The 



reason for this is the following. To a good approximation a composite of k nucleons is an 
incompressible sphere with volume /c/po where the value of po is — 0.16 fm"^. The volume 
available for translational motion (eq.(3)) is then V = V — Vexduded where we approximate 
Vexciuded — N/ Pq — Vq the normal volume of a nucleus with N nucleons. Similar corrections 
are implicit in Van der Waals equation of state. This is meant to take care of hard sphere 
interactions between different particles. This answer is approximate. The correct answer is 
multiplicity dependent. The approximation of non-interacting composites in a volume gets 
to be worse as the volume decreases. We restrict our calculation to volumes V greater than 
2Vo. This is how the calculations reported in the next section proceed. We choose a value 
of Vq/V — pI Pq from which V^jV — pj p^ — p/{po — p) is deduced. This value of p is used 
in eq.(8) to calculate p and all other quantities. We plot results as function of p/po- If we 
plotted them as function of p/po the plot would shift to the right. 

Calculations in [3] are continued beyond the limit p/ po = 0.5. They find one can identify a 
critical point at T = Tc = IS.OMeF, p/po — 1 and Pc — oo. At very high pressure the model 
should break down: Zk{int) must change at such high pressure although nuclear physics 
says that nuclei being highly incompressible moderate pressures should leave the internal 
partition functions relatively unchanged. Another way of saying this is that interaction 
between composites should be taken into account for ~ Vq. 

III. THE P-p CURVES FOR ISOTHERMALS 

For a given temperature T and p we solve for p and then pressure. This is plotted in 
Fig. 1. For each isotherm shown the pressure rises rapidly at first with p but then flattens 
out. The flattening depends upon the value of kmax- For low kmax (shown in section VII) 
there will not be any flattening. The value of k^ax used in Fig. 1 is 2000. There is still 
a very slight rise in p (not discernible in the Figure). The figure empirically allows us to 
designate two regions: a purely gas phase where the pressure rises with density and a liquid- 
gas co-existence phase where the volume changes but the pressure is nearly stationary. One 
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way of seeing this is that p — T{M/V). As V decreases so does M so as to compensate in 
the co-existence region. Fig. 1 also shows that in the gas phase, the chemical potential rises 
rapidly with density but then flattens out in the co-existence phase. 

A discussion about /i dependence in the co-existence region is in order here. For k^ax 
large p is much more sensitive to than pressure (p has weighting of fc^^^, eq.(8) whereas p 
has A;^/^ weighting (eqs. (5) and (7)). For very large kmax an inflnitesimal change in /j, will 
lead to a finite change in p but only a very small change in p. In the limit k^ax — > oo we 
will reach ideal liquid gas phase transition: no change of p in the coexistence region and no 
change of pressure. This is demonstrated in section VI. 

IV. WHAT CONSTITUTES THE GAS AND WHAT CONSTITUTES THE 

LIQUID? 

As an example, at fixed temperature 7 MeV we show in Fig. 2 the distribution of compos- 
ites (a) in the pure gas region {p/po — 0.12)) and (b) in the co-existence region {p/po — 0.22). 
In the gas phase the sum rule J2knk/N = 1 is exhausted well before we reach 50. There 
are no heavy composites. In the co-existence region there are light particles {k < 40), then 
nothing for a long range of k, and then there are heavy particles with k between 1800 to 
kmax = 2000 (the figure shows the population of fc=1900 to k^ax = 2000). A safe functional 
definition for the gas phase is all composites between k=l to k = 100 and for the liquid 
phase all composites between k^ax ~ 300 to kmax- Thus both the liquid and the gas phases 
are quite complicated; consisting of not one or two but many species although they are all 
made up of the same elemental nucleon. 

V. THE CLAUSIUS-CLAPEYRON RELATIONS 

In the co-existence region the pressure is a rapidly increasing function of the temperature. 
In Fig. 1 these are the fiat regions shown for three temperatures (dash-dot 6 MeV, solid 7 
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MeV and dash 7.5 MeV). The Clausius-Clapeyron relation for hquid-gas phase transition 
provides an equation for the rate of change [7]: 

dT Av 

where As can be taken to be the change of entropy per unit mass and Av the corresponding 
change in volume as matter moves across phase transition. We can take these changes to 
be per nucleon. The following substitutions are made:As = L/T (L is the latent heat) and 
Av — Vgas — viiq. The standard approximation now is Vgas » vuq^Vgas — ^/ Pgas — T/p and 
thus 

(10) 

If we make the assumption that L is nearly independent of temperature then the equation 
integrates out to give 

Inp^lnpo- L/T (11) 

This does not work well in our case (Fig. 3): In p is not a hnear function ol (3 — 1/T. We 
can discard the assumption that L is constant and instead use eq.(lO) to get an idea of L 
using values of dp/dT,p and T from Fig. 1. If this is done then at 6 MeV temperature, the 
value of L turns out to be 54 MeV and at 7 MeV temperature this reaches 70 MeV. Since 
the binding energy per particle for an infinite cluster is 16 MeV, these values are clearly 
unacceptably high. 

Let us ask what went wrong in going from eq.(9) to eq.(ll). The passage from eq.(lO) to 
eq.(ll) assumed that the latent heat is independent of temperature. We will show that this 
is approximately correct. However, the approximation Vgas — T/p is very inaccurate and 
depending upon the temperature, corrections due to vug can be significant. When all this is 
taken into account eq.(9) is satisfied remarkably well. We elaborate first on the latent heat. 
For definiteness fix on the isothermal at 7 MeV. From Fig. 1 we can determine the density at 
which the system enters the purely gaseous phase and its energy per particle from J2 nkCk/N 
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where = (3/2)T + ek{int). Here ek{int) consists of volume energy (which is negative), 
surface energy and as well contributions from excited states. The expressions are given in 
section II. The nucleons are passing from a liquid state (from a composite with k ~ 2000) to 
the gaseous phase. The energy per particle in the liquid phase was taken from a composite 
of 1950 particles. Reasonable variation around this number will only change the calculated 
value slightly. Latent heat per particle calculated is 12.66 MeV at temperature 6 MeV and 
11.55 MeV at 7.5 MeV. 

By far the major error is in assuming that Vgas — T/p. The pressure is given by p = 
{M/V)T and not {N/V)T where M is the multiplicity and of course the total number 
of particles. Thus p — {M/N)T/vgas- The factor (M/N) when the system just turns into a 
pure gas phase is 0.276 at T=6.0 Mev, 0.194 at 7.0 MeV and 0.152 at 7.5 MeV. Writing a 
for M/N we find that eq.(lO) should be rewritten as 

^ = ^ (12) 

where we have used the fact that vuq = 1/po- As an example, at 7 MeV ^ = 0.0535fm~^ 
from Fig. 1 and 0.0530fm~^ estimated from formula (12). At 6 MeV the corresponding 
numbers are 0.0129fm~^ and 0.0108fm~^ respectively. 

VI. THE LIMIT LARGE Kmax- 

We will now consider changes in the values of various quantities as we change from one 
large value of kmax to another large value of kmax- For definiteness we will concentrate on 
one isothermal, for example the T=7 MeV case. In Fig. 1 we have two regions: pure gas 
phase and the co-existence phase. The pure gas phase is trivial. Nothing changes as we go 
from one large kmax to another; kmax— 2000 is large enough in this case. It is easy to see 
why results become insensitive to changes in the value of kmax- In the gas phase there is no 
population in high k composites so it does not matter whether the summation stops at a 
given high value of kmax or another high value of kmax- The situation is more complicated 
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but also more interesting in the co-existence region as we have population both at the lower 
end and the higher end of k. 

In Fig. 4 we plot the values of jj, and pressure calculated for kmax in the range k^nax — 2000 
to 5000. This is done at a fixed value of p/po — 0.3 which is in the co-existence region. As 
shown in the figure, both curves are well fit by a parametrisation a + bexp{—ckmax) where 
kmax > 2000, with values given in the caption. This means that within the accuracy with 
which this calculation was carried out, the values of p, and pressure in the infinite kmax limit 
are -18.504 MeV and 0.0294 MeV fm"^ respectively 

A similar calculation as above was done for p/po=0.4. The fitted values of a, b and c for 
p were -18.504 MeV, 0.33892 MeV and 0.0003197 respectively. For pressure, the parameters 
were 0.0294 MeV fm-^ 0.007633 MeV fm-^ and 0.0003906 respectively. Note that the 
extrapolation demonstrates that that neither p nor the pressure change in the co-existence 
region in the limit kmax oo. This firmly establishes the present model as a model of 
liquid-gas phase transition as was stated in section III. 

Lastly we want to estabish a scaling law. Given the fractional occupation fk{kmax) = 
kuk/N for a large value of kmax, do we know the fractional occupation fk'{k'max) another 
large value of kmax'^ Based on the discussion, so far we expect that if one is in the purely 
gas phase fkikmax) = fk{k'max) this is indeed the case. 

In the co-existence phase a lowest order approximation is based on the following ap- 
proximation. We expect the fractional occupation to match near the beginning {k small), 
near the end (near k < k-^ax and k' < k'^^^) and in between there is almost no occupation. 
Thus for k small fk{kmax) ^ fk{k'max) and near the high end fk{kmax) ^ fk'{k'max) where 
kmax — k — k'^^^ — k' . This is not very accurate but an accurate representation for low k is 
given using the parametrisation: 

k 

In fkikmax) = In fkik'^ax) + ^b[exp{-ckmax) - exp{-ck'^^^)]. (13) 

An equation relating the large clusters can also be written down, but the functional form is 
quite complicated. 
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VII. SMALL SYSTEMS: AN EXACT CANONICAL MODEL SOLUTION 



The model can be solved when the number of particles N is finite. Extensive use of 
the canonical model has been made to fit experimental data [8] so just an outline will be 
presented for completeness. Among other applications, the canonical model can be used to 
study finite particle number effects on phase transition characteristics. 

Consider again TV identical particles in an enclosure V and temperature T. These 
nucleons will combine into monomers, dimers, trimers etc. The partition function of the 
system in the canonical ensemble can be written as 

Q-=En^ (14) 

Here Zi is the one particle partition function of a composite which has i nucleons. We already 
encountered Zi in section II: Zi — Zi{tran)zi{int) with Zi{tran) and Zi{int) given in detail. 
Other forms for Zi can be used in the method outlined here. The summation in eq.(14) is 
over all partitions which satisfy A" = X^irij. The summation is non-trivial as the number of 
partitions which satisfy the sum is enormous. We can define a given allowed partition to be 
a channel. The probablity of the occurrence of a given channel P{n) = P(ni, n2, ns....) is 

m = ^W-^- (15) 

Wn nil 

The average number of composites of i nucleons is easily seen from the above equation to 
be 

(rii) = z~-=^ (16) 

VAT 

Since Yli^'^i = one readily arrives at a recursion relation [9] 

1 ^ 

Qn — kzkQN-k (17) 

k=l 

For one kind of particle, Qn above is easily evaluated on a computer for A" as large as 
3000 in matter of seconds. It is this recursion relation that makes the computation so easy 
in the model. Of course, once one has the partition function all relevant thermodynamic 
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quantities can be computed. For example, eq. (7) still gives the expression for pressure 
although one could correct for the center of mass motion by reducing the multiplicity by 1: 
p = T(M — 1)/V . The chemical potential can be calculated from fi = F{N) — F{N — 1) 
where the free energy is F{N) — —T In Qn which is readily available from the calculation. 

In Fig. 5 we show an example of the canonical model calculation. The temperature is 
6 MeV. The number of particles N is 200. The value of the largest allowed cluster kmax is 
also 200. This would be a typical nuclear physics case. In the same figure we also show the 
result of a grand canonical calculation with the same kmax (of course for the grand canonical 
N is very large). At small density the results are the same but they become different at 
larger densities. In the canonical result there is a small region where dp/ dp is negative. 
This is a finite particle number effect since for large (grand canonical result) any negative 
compressibilty disappears. Negative compressibility can lead to negative value for Cp [8]. 
The grand canonical result shows that for kmax—'^^O typical liquid-gas co-existence is not 
found and there is no region where p is constant when the density changes. 

VIII. DISCUSSION 

Results in section III to section VI show that the model of excited matter breaking up 
into clusters with saturation properties leads to a first order phase transition. This has 
relevance to heavy ion collisions at intermediate energy but may have significance in other 
areas of physics as well. This model for first order phase transitions is extremely easy to 
implement. A very significant advantage of the model is that it can be solved not only in 
the thermodynamic limit (large A^) but also for a finite number of particles. Thus one can 
study how observables change as one progresses from small to large systems. 

We like to end this discussion by noting that in spite of a very different approach that 
is adopted here to arrive at the key equations 5, 7 and 8, formally the pressure and density 
equations have the same structure as those encountered in the well-known Mayer cluster ex- 
pansion [10]. These are p = ^^^"^P'^' E^i exp{(3p,k)hk and p = (HH^p!^! y:T=i exp{Pfik)kbk. 
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Here instead of the cluster integral we have k^/'^Zk{int). 
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FIG. 1. Behavior of pressure p and chemical potential ii against p//9o for 3 different temper- 
atures : dash (7.5 MeV), solid (7 MeV) and dash-dot (6 MeV). We identify as purely gas phase 
the region where the pressure and the chemical potential [i rise with density and the co-existence 
region where they remain constsnt. 
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FIG. 2. At temparature 7 MeV yields of composites at two densities p/po=0-12 (gas phase) and 
p/Po=0.22 (the co-existence phase). For the first case there are no heavy composites (dashed hne). 
The sum rule ^kuk/N = 1 is already satisfied to good accuracy by A; = 40. In the second case 
(solid line) there are light particles (less than 50 nucleons) and there are heavy particles (greater 
than 1800 particles). Together these exhaust the sum rule. In k space there is a huge gap for 
particles between large and small. The occupation number in this region is very close to 0. 
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FIG. 3. A plot of ln{p) against the inverse of temperature. The relationship is not linear. Here 
Po is 1 MeV fm~^. 
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FIG. 4. The solid curve in the upper panel is a plot of jj. against kmax in the range of kmax='^^0^ 
to 5000 with density at p/po=0.3 and temperature 7 McV. The dashed curve is a fit with the 
parametrisation a + bexp{—ckmax)- The values of the fit parameters are a = —18.504 MeV, 
6=0.33748 MeV and c=0. 0003842. Similar quantities for pressure in the lower panel. The fit 
parameters are a=0.0294 MeV fm^^, 5=0.007503 McV fm"''^ and c=0.0003842. Similar curves for 
p/po=0.4 yield equally good fit and give the same values for a but different values for b and c. 
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FIG. 5. The solid curve in the upper panel is a plot of pressure against density in the canonical 
model. The number of particles is A^=200 exactly and kmax is also 200. Note that there is a region 
of negative compressbilty. The dashed curve is the grand canonical result with the same k^ax ■ The 
two curves coincide at low density. Note that in the grand canonical model the increase of pressure 
with density goes down later but never disappears for this low kmax- The lower panel compares 
the chemical potentials. 
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